# konrad.smolinski@gmail.com
# date: 		26/11/2010
# last update: 	12/05/2011
#
# info:
#	this file works with an illustration from rcprobit4.pdf and rcprobit1.nb

# 	following 'rcprobit4.pdf'

#	(adc) mail note: 13/05/2011
#	In my initial grid experiments I have used a binary instrument 
#	z in {-1,1}  
#	a1=0, a2=1, b1=1, b2=1  and d1=1 (not 2).  
# ------------------------------------------------------------------
# Call for libraries:
rm(list=ls())
source("./lvlLibSrc.R")
source("./rcLibSrc.R")

# CONSTRUCTOR - search (a0,a1) parameters
rcCrFunDef <- function(  a 		= c(0,1) 
						,b 		= c(0,1)
						,bpar	= c(0,1)
						,d 		= c(1,3^(-1/2),3^(-1/2),3^(-1/2))
						,zval 	= c(-0.5,0.5)
						,xthres	= c(0)
						,xval	= c(0,1)
						,acc 	= 40
						){
							
	K <- length(xval)
	R <- length(zval)
	dgp <- rcDgp(a,b,d,xthres,xval,zval)
	
	# this part is done once:
	probIdx <- rcLeaveOutToProbs(K)
	probMat <- matrix(0,ncol=R,nrow=length(probIdx))

	for(r in 1:R){
		for(i in 1:length(probIdx)){
			 probMat[i,r] <- sum(dgp[[r]][probIdx[[i]]]) 
		}
	}
	rhs <- apply(probMat,1,max)	


	
	ff <- function(par){
		xval <- c(0,1)		
		K <- length(xval)		
		# integrand over Base Components - this part depends on initial parameters
		loIdx <- rcLeaveOutIdx(K)
		lhs <- rep(0,length(loIdx))
		
		# parameters enter here:
		intVec <- rcBaseIntegrate(a=par,b=bpar,xval=xval)
		for(i in 1:length(loIdx))	lhs[i] <- sum(intVec[loIdx[[i]]])

#		lhs > rhs
		return(all(lhs > rhs))		
	}
	
return(ff)	
}

# ---------
b1 <- c(-0.4,-0.2,0,0.2,0.4,0.6)

resList <- list()
# Search
for(i in length(b1)){
	rcff <- rcCrFunDef(b=c(b1[i],1))
#	plot(c(-1,1),c(-1,2),xlab=expression(a[1]),ylab=expression(a[2]),col="white")
	res <- lvlFind2D(c(0,1),criterion=rcff,step=c(0.02,0.02),rderr=0.01,plt=F)
#	points(res[1],res[2],col=2,pch=19,cex=0.7)
#	polygon(res[,1],res[,2],col=2)
	resList[[i]] <- res	
}


save(resList,file='../Data/rcEx6-2-Data.RData')
load(file='../Data/rcEx6-2-Data.RData')


# PLOT 
sv <- T

for(i in 1:length(res)){
	res <- resList[[i]]
	xlattice <- unique(res[,1])
	ylattice <- unique(res[,2])

	xrange <- diff(range(res[,1]))
	yrange <- diff(range(res[,2]))

	xmin <- min(res[,1]) - 0.1*xrange
	xmax <- max(res[,1]) + 0.1*xrange
	ymin <- min(res[,2]) - 0.1*yrange
	ymax <- max(res[,2]) + 0.1*yrange

	if(sv) pdf(paste("../Graphs/rcEx6-2-",i,".pdf",sep=""),width=10,height=6)
		plot(c(xmin,xmax),c(ymin,ymax),xlab=expression(a[1]),ylab=expression(a[2]),col="white")
		polygon(res[,1],res[,2],col="yellow",border=1)
		points(res[,1],res[,2],col=1,pch=19,cex=0.3)
		points(0,1,col=2,pch=19,cx=1.2)
		abline(h=1,col=2)
		abline(v=0,col=2)
	
		text(xmax-0.2*xrange,ymax-0.1*yrange,expression(b[1]==0))
		text(xmax-0.2*xrange,ymax-0.2*yrange,expression(b[2]==1))
	if(sv) dev.off()

}